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Convection  Algorithms 
for  Magnetospheric  Particles 


1,  INTRODUCTION 

This  report  documents  in  detail  the  implementation  of  the  convection  model  used 
in  the  research  conducted  by  Kerns,  Hardy,  and  Gussenhoven  [1993],  The  convection 
model  assumes  an  Earth-centered  dipole  magnetic  field  with  equipotential  field  lines,  and 
it  assumes  a  static  cross  tail  electric  field  of  the  form  VoeL^sin(<^)  where  y  is  a  positive 
real  number,  L  is  L-shell  fsee  Mcllwain  [1961]),  and  ^  is  local  time.  Convection  occurs 
because  of  ExB  drifts  and  magnetic  field  gradient  and  curvature  drifts.  The  model  is 
used  to  calculate  drift  paths  of  particles  as  well  as  to  calculate  the  range  of  kinetic  energy 
for  which  ions  and  electrons  move  on  drift  paths  open  to  the  magnetospheric  tail.  Drift 
paths  are  calculated  using  conservation  of  the  first  two  adiabatic  invariants  [Chen,  1984] 
and  conservation  of  energy.  Electrons  are  shown  to  be  on  open  drift  paths  if  the  kinetic 
energy,  W,  is  0<  W<  where  is  the  electron  cutoff  energy,  and  ions  are  on  open 
drift  paths  if  W^<  W<  Wj,  where  and  are  the  ion  lower  and  upper  cutoff  energies 
respectively.  The  algorithms  presented  in  this  document  are  optimized  to  converge  to  a 
solution  quickly  for  any  positive  real  value  of  y,  for  any  particle  energy,  position,  and 
equatorial  pitch  angle,  and  for  any  cross-tail  field  strength.  These  algorithms  are  used 
by  the  CRRES  survey  software  called  "SHOWLEPA.EXE"  which  was  developed  by  this 
author  at  the  Phillips  Laboratory  in  support  of  the  research  by  Kerns,  Hardy,  and 
Gussenhoven  [1993]. 

Implementations  of  equations  are  given  in  algorithms  that  are  listed  in  text  boxes 
separate  from  the  main  text.  These  algorithms  are  written  in  pseudo-code  similar  to 
Pascal  to  make  them  readable  to  the  non-programmer  and  useable  to  the  programmer; 
comments  are  given  in  italics  and  are  not  considered  part  of  the  instructions  or  code. 
These  are  added  to  make  the  code  more  readable.  Constants  are  defined  at  the  top  of  any 
procedure  or  function  using  the  "  = "  symbol,  and  each  constant  represents  a  value  that 
does  not  change.  Constants  are  used  to  reduce  unnecessary  math  to  make  the  code 
execute  more  quickly.  All  variables  are  assumed  to  contain  real  values  unless  specifically 
identified  as  integers,  and  the  variables  required  for  any  given  function  or  procedure  are 
listed  near  the  top  of  the  code  for  that  procedure  or  function.  Algorithms  use  logical,  not 
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mathematical,  statements  so  that  the  statement  "i=i+r',  which  is  nonsense  in 
mathematics,  means  to  increment  the  vaiue  contained  in  the  variable  called  "i”  by  1. 
Subscripts  are  part  of  variable  names,  and  are  used  to  make  code  more  readable.  TTie 
one  exception  to  the  use  of  subscripts  is  with  the  table  of  values  for  the  function  I(y) 
given  in  Algorithm  1;  here  subscripts  are  used  to  identify  individual  elements  of  the 
table.  Function  and  procedure  names  are  given  in  bold  letters,  and  all  parameters  passed 
to  these  use  the  naming  convention  given  in  Algorithm  1.  The  repeat  block  until 
condition  loop  implements  a  block  of  code  once,  and  then  continues  to  re-implement  the 
same  block  until  the  given  condition  is  true.  The  while  condition  block  iterates  the  next 
block  of  code  for  as  long  as  the  given  condition  is  true.  The  if  condition  then  code 
statement  implements  the  block  of  code  only  if  the  condition  is  true,  and  the  if  condition 
then  codel  else  code2  statement  implements  the  "code!"  block  of  code  if  the  condition 
is  true  and  the  "code2"  block  of  code  otherwise.  All  blocks  of  code  that  are  longer  than 
one  line,  with  the  exception  of  those  bounded  by  "repeat  .  .  .  until"  loops,  are  marked 
by  "begin"  and  "end"  statements  .  Indentation  is  used  to  make  loops  and  blocks  more 
readable.  Nested  parentheses  are  made  more  readable  by  using  "{[(..)]}". 


2.  CONVECTION  MODEL 


In  this  section  the  model  used  to  calculate  the  trajectories  of  central  plasma  sheet 

Algorithm  1  Global  constants  and  variables  used  by  other  algorithms,  and  names  of  parameters 
passed  to  functions  and  procedures  given  in  other  algorithms. 


Global  constants  used  in  all  calculations  and  algorithms; 


uj  s  7.272X10® 
a  =  6.378X10* 

Bo  s  3.1X10® 

I,  s  l(i/N);  O^i^N 
N  s  1000 
6^  =  0.01 
tfw  =  0.001 
=  10000 


radians  per  second.  Earth 's  spin  rate 

meters.  Earth’s  radius 

tesia,  magnetic  field  at  Earth's  Equator 

to  defined  in  Algorithm  2;  "/"  is  element  index  of  table. 

is  first  element,  and  1^  is  last  element  of  table. 

Rg,  error  allowed  in  calculating  L-Shells 

eV,  error  allowed  in  calculating  cutoff  energies 

eV,  energy  step  used  in  cutoff  solutions 


Global  variables  used  in  ail  calculations  and  algoritlims; 
A  volts/meter,  from  Equation  2 

Y  from  Equation  2 

Lo  from  Equation  7 


Standardized 

y 

s 

W 

L 

d 

J 

U 


names  of  parameters  passed  to  functions  and  procedures: 
sine  of  equatorial  pitch  angle,  (0.0  ^  y  s  1.0),  (Equations  3  and  4) 
sine  of  local  time;  sin(n/12-LT),  (-1.0  s  s  <  1.0) 
kinetic  energy  in  eV 
L -shell  in  Rg 

first  adiabatic  invariant  constant  (Equation  3) 
second  adiabatic  invariant  constant  (Equation  4) 
total  energy  of  particle  in  eV  (Equation  5) 


electron  boolean;  true  if  particle  is  an  electron  and  false  if  it  is  an  ion. 
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(CPS)  particles  in  the  inner  magnetosphere  under  conditions  of  steady  convection  is 
described.  The  model  follows  the  same  method  and  assumptions  used  by  Volland  [1973]. 
It  assumes  a  dipole  geomagnetic  field  aligned  with  the  Earth’s  spin  axis  and  having 
equipotential  field  lines.  The  co-rotation  electric  field  potential  is  then  taken  as 

V  =  Ba^o3L'\  (1) 

itt  0  ^ 

where  Bg  is  the  strength  of  the  equatorial  magnetic  field  at  the  surface  of  the  Earth,  a  is 
the  Earth’s  radius,  u  is  the  Earth’s  spin  rate,  and  L  is  L-shell  in  Rg.  The  cross-tail 
potential  is  assumed  to  be  of  the  form  given  by  Volland  [1973] 

V  =  aAL^s\niil>)  (2) 

where  A  scales  the  field  strength,  y  determines  the  rate  of  decrease  in  potential  with 
decreasing  altitudes,  and  0  is  the  local  time.  The  model  is  used  to  determine  the  range 
of  kinetic  energy  of  ions  and  electrons  in  which  .the  particles  convect  on  open  drift -paths 
from  the  tail.  The  bounds  of  this  range  are  referred  to  as  cutoff  energies,  and  these  are 
a  function  of  the  position  at  which  the  ^article  is  measured  as  well  as  the  parameters  A 
and  y  of  Eq.  (2).  Cutoff  energies  for  steady  convection  have  been  introduced  previously 
in  other  works  including  Stem  [1975]  and  Southwood  and  Kaye  [1979].  In  all  examples 
given  in  this  report,  7=2,  after  Stem  [1975],  and  the  values  of  the  constants  in  Eqs.  (1) 
and  (2)  are  given  in  Algorithm  1 .  The  values  of  A  and  y  are  given  as  global  variables 
in  Algorithm  1,  and  these  are  used  by  the  subsequent  algorithms.  These  values  are 
scaled  to  represent  different  magnetospheric  conditions.  In  addition  to  the  values  of  the 
constants  used.  Algorithm  1  contains  the  standard  names  of  parameters  passed  to  the 
functions  given  in  subsequent  algorithms. 


Algorithm  2  Implementation  of  I(y)  given  in  Equation  4. 


FUNCTION  DEFINITION  OF  My) 
local  constants  required  for  Uy): 
a  =  1.380173 

s  0.055  32  *  -0.037 

bo  3  2a  b,  s  3a, 

b4  s  -6a4  b^  s  2(8-83) 
local  variables  required  for  l(y); 

Yl»  Vc»  1 

implementation  of  l(y): 
if  y  =  0.0  then 
I  =  bo 

else  if  y  =  1.0  then 
I  =  0.0 


S  s  y,a-  0.3702402 

83  s  -0.074  84  s  0.056 

bj  s  682  bj  ®  4B-2a-3a,-6a2-t-6a4 

b,t  SE  -48 


else 

begin  conditional  block 

Yl  =  Inly)  Inixl  is  natural  logarithm  of  X 

Yc  ~  expOAy^)  expfx}  is  inverse  natural  logarithm  or  ef 

I  =  bo  +  {b,  +  (bj  +  (bg  +  bi-yL+  b4’yc)-yc]-yc}-Vc  +  b^-y’^ 
end  conditional  block 
return  value  of  1 
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Particle  motion  is  assumed  to  be  a  combination  of  ExB  drift  and  magnetic  field 
gradient  and  curvature  drift.  Trajectories  are  determined  by  conservation  of  energy  and 
conservation  of  the  first  two  adiabatic  invariants.  From  these  the  three  following 
constants  of  motion  are  defined: 

H  =  y^WV  (3) 

J  =  W'LKy) 

U  =  W q[aAL^sini<t>)  -L-%a“u)]  (5) 

where  n  and  J  are  proportional  to  the  first  two  adiabatic  invariants,  U  is  the  total  energy, 
y  is  the  sine  of  the  equatorial  pitch  angle,  W  is  kinetic  energy,  and  q  is  charge.  Ions  are 
assumed  to  be  singly  ionized  so  that  q=±i  for  ions  and  electrons  respectively.  The 
function  I(y)  in  Eq.  (4)  is  an  integration  function  of  the  bounce  motion,  and  the 
approximation  given  by  Ejiri  [1978]  is  used  for  the  implementation  shown  in  Algorithm 
2.  Equation  (5)  is  the  sum  of  the  kinetic  energy  and  potential  energy  [Eq.  (1)  and  (2)]. 


3.  CAI  CULATION  OF  DRIFT  PATHS 

Valid  drift  paths  of  particles  are  the  set  of  all  points  where  Eqs.  (3),  (4),  and  (5) 
give  and  U=U^f.  These  constants  of  motion  are  for  a  particle  with 

measured  kinetic  energy  W=W^f  at  a  position  and  y=yv. 

Algorithm  3  shows  the  implementation  of  these  equations.  Here  the  function  I(y)  from 
Eq.  (4)  is  replaced  by  a  look-up  table  with  N+1  equally  spaced  elements,  and  the  value 
of  Ify)  is  found  by  interpolating  between  the  two  nearest  table  elements.  This  table  is 
used  to  make  the  algorithms  faster. 

The  simplest  case  of  convection  is  for  zero  kinetic  energy  particles  (Wj^=0),  and 
this  is  presented  first.  From  Eqs.  (3)  and  (4),  and  conservation  of  energy 

[Eq.  (5)]  requires  that  these  particles  drift  on  equipotential  lines  given  by 

V  =  aAL^sm{<fi)-L'^Bga^(>)  (<S) 


Figure  1  shows  the  equipotentials  given  by  this  equation  for  /4=7.46xl0’*  V/m.  In  this 
figure,  noon  is  towards  the  top  of  the  page,  and  each  hour  of  LT  is  marked  by  a  radial 
line.  The  half  shaded  circle  represents  the  Earth,  and  the  concentric  circles  show  L-shells 
2  through  9  Re-  The  heavy  lines  show  equipotentials  from  -45  kV  to  5  kV  at  5  kV 
intervals,  and  these  are  symmetric  about  the  dawn-dusk  line.  In  this  example,  the  three 
lines  shown  closest  to  the  Farth  (-45  to  -35  kV)  are  closed,  and  zero  kinetic  energy 
particles  on  these  lines  are  trapped  and  drift  around  the  Earth  in  an  Eastward  direction. 
The  equipotential  lines  further  from  the  Earth  are  open,  and  zero  kinetic  energy  particles 
on  these  are  not  trapped,  but  drift  sunward  from  the  tail  until  they  encounter  the  mag¬ 
netopause. 
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A^orithm  3  Implementation  of  Equations  3,  4,  and  5. 


FUNCTION  DEFINITION  OF  p(y,  W.  LI  Equation  3. 

implementation  of //(y,  W,  L): 

return  value  of  |yL)*-L-W 

FUNCTION  DEFINITION  OF  J(y,  W.  L)  Equation  4. 

local  variables  required  for  J{y,  W,  L): 

I,  value  of  KyJ  interpolated  from  the  array  /,  (Algorithm  1) 

i  index  of  element  in  i,  to  be  used  for  interpolating  /, 

implementation  of  JW-  W,  L); 
if  y  =  0.0  then 
I,  =  Iq 

else  if  y  =  1.0  then 

I,  =  In 

else 

begin  conditional  block 

I,  =  N  y  N  is  maximum  index  of !  (Algorithm  1) 

i  =  truncd,)  trunc  returns  the  largest  integer  less  than 

I,  =  t,  -  i  or  equal  to  /, 

I,  =  l|‘(1-l,)  +  U,  l,  now  contains  interpolated  value  of  Ky) 
end  conditional  block 
return  value  of  W’‘'L  I, 


FUNCTION  DEFINITION  OF  Ulelectron,  W,  L.  s) 
local  variables  required  for  Uielectron,  W,  L,  s); 

k  5E  Bo-a*-w 

implementation  of  Ulelectron,  W,  L,  s): 
if  electron  then 

return  value  of  W  -  a-A-s-U'  +  k/L 

else 

return  value  of  W  +  a-A  s-L’'  -  k/L 


Equation  5. 


see  Algorithm  1 
q  =  -1 

0=  +  1 


The  zero  energy  Alfv^n  layer,  is  the  last  closed  drift  path  for  zero  kinetic 
energy  particles.  In  Figure  I,  tear-drop  shaped  and  lies  along  the  -30  kV 
equipotential,  and  it  bounds  the  shaded  region  Aat  surrounds  the  Earth.  at  dusk 
(1800)  is  a  stagnation  point  for  W=0  particles  because  here  the  corotation  electric  field 
is  equal  and  opposite  to  the  cross-tail  field  resulting  in  no  electric  field  drift;  Ejiri  [1978] 
derives  the  L-shell  of  this  point,  1^,  as: 


aaB 

o 

7^ 


(7) 


The  value  of  Z,  is  used,  instead  of  the  value  of  .4,  to  indicate  the  cross-tail  field  strength 
in  the  remaining  examples  because  .d  can  be  determined  from  Eq.  (7),  and  its  value  scales 
with  the  distance  of  %  from  the  Earth.  Zo=4.58  Rg  in  Figure  I . 

The  L-shell  of  can  be  determined  for  any  local  time  by  using  Eqs.  (6)  and  (7). 
The  potential,  V^,  of  %  is  found  by  substituting  4)  for  L,  and  -1  for  sin(<^)  in  Eq.  (6). 
As  is  evident  in  Figure  I,  is  at  maximum  L-shell  at  dusk,  so  that  Lq  is  always  greater 
than  or  equal  to  the  L-shell  of  Xo  for  the  given  local  time.  Solving  for  the  partial  of  V 
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00 


Figure  1  Equipotentials  for  i4=7. 46x10’  V/m  [Eq.  (2)].  Radial  distance  shows  L-shell  in  Rg, 
and  angular  distance  shows  local  time.  The  shaded  region  is  bounded  by 


[Eq.  (6)]  with  respect  to  L  and  setting  the  result  equal  to  zero  gives 


■^av/ar-o 


—yAsmi4>) 


I 


Y^l 


(8) 


Tnis  equation  shows  that  the  slope  of  V  as  a  function  of  L  in  Eq.  (7)  does  not  change  sign 
on  the  dawn  side  (sin(<f))^0)  because  dV/dL  cannot  equal  zero.  On  the  dusk  side 
(sin(^)<0),  Eq.  (8)  shows  that  the  L-shell  at  which  the  slope  changes  sign  (dVldL=0) 
must  be  greater  than  [Eq.  (7)].  Because  the  slope  does  not  change  sign  in  the  region 
where  the  valid  solution  exists  O^L^Lo),  a  simple  binary  search  can  be  used  to 
determine  the  L-shell  at  a  given  local  time  with  the  same  potential  as  that  for  at  dusk. 
The  starting  limits  of  L-shell  for  the  search  are  0  and  1^,  and  the  implementation  is 
shown  in  Algorithm  4.  This  algorithm  assumes  that  the  values  of  A,  y,  and  have  been 
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Algorithm  4  Calculation  of  the  L-shell  of  for  a  local  time  specified  by  sin(^). 


FUNCTION  DEFINITION  OF  L*{s) 
toc3^  nstants  required  for  L^ls): 

k  = 

local  variables  required  for  L^ls): 

k,,  V,  V<{/  ALg,  Lg 
implementation  of  Lg(s): 

Vg  =  k/Lo  -  a-A  Lo'' 
k,  =  a-A-s 
Lg  =  Lo 
ALg  =  Lo 

if  s> -0.999999  then 


potential  of  «  found  from  Equation  6 
These  multiplications  are  done  before  the  loop 
highest  possible  value  for  ig 

step  size  assuming  lowest  possible  value  of  ig  is  zero. 
Prevent  round  off  error  from  causing  a  run-away  loop. 


repeat  the  following  loop: 


V  =  k/U  + 


ALg  = 
if  IV  > 

else 


k.Lg-' 


’/^■ALg 

Vg)  then 

l-a  ~  Lg  -  ALg 


potential  from  Equation  6  at  sin(0J  =  s 
the  step  size  halves  during  each  iteration 

Ig  should  be  smaller 


until  ALg  < 
return  value  of  Lg 


Lg  =  Lg  +  ALg  fg  should  be  larger 

The  hop  is  terminated  once  sufficient  accuracy 
is  reached,  and  the  value  of  is  returned. 


previously  determined.  A  similar  algorithm  not  shown  in  this  report  was  used  to  plot  the 
equipotential  lines  in  Figure  1.  Algorithm  5  shows  the  implementation  of  Eq.  (7)  to  set 
the  value  of  the  global  variable  L^,  when  A  and  y  are  known.  It  also  shows  a  simple 
method  of  determining  the  values  of  A)  and  A  when  only  y  is  known  so  that  includes 
a  specified  L-shell  and  local  time.  It  uses  a  binary  search  to  determine  the  value  of  A 
that  results  in  a  value  of  Lg(s)  (Algorithm  4)  equal  to  L. 

When  W>0,  Eqs.  (3),  (4),  and  (5)  must  be  solved  simultaneously  to  determine 
the  drift  paths  of  a  particle  measured  with  the  constants  of  motion,  n^,  Ju,  and  Um. 
Equations  (3)  and  (4)  are  first  combined  to  solve  for  y  as  a  function  of  L-shell. 

l(y)y-^ =  Q  (9) 

Equation  (9)  can  be  solved  analytically  if  the  function  I(y)  is  assumed  to  be  linear  such 
that  I(y)  »  a,  -I-  by.  The  solution  of  y  is  then 

y  *  a,iQ-b,r 


The  values  of  a,  and  b,  are  determined  using  the  two  elements  of  table  I  (Algorithm  1) 
nearest  to  the  solution  of  y.  This  method  leads  to  an  iterative  solution  which  is  illustrated 
in  Figure  2.  The  two  heavy  lines  show  the  function  I(y)  and  yQ,  where  Q=l.69.  The 
intersection  of  these  two  lines  is  the  solution  to  Eq.  (9).  At  the  beginning  of  the 
iteration,  y  is  assumed  to  equal  1.0,  and  I(y)  is  approximated  by  a  line  at  Point  A.  This 
approximation  is  shown  by  the  thin  line  through  Point  A.  Equation  (10)  gives  y=0.467 
at  the  intersection  of  this  ^proximation  with  yQ.  A  new  line  is  used  to  approximate  I(y) 
at  Point  B  using  the  new  v^ue  for  y.  This  approximation  yields  Point  C  at  y =0.5 12,  and 
then  converges  on  the  next  iteration  to  Point  D  at  y =0.5 13.  W  is  then  determined  using 
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Aigorithm  5  Calculation  of  the  global  variables  and  A. 

PROCEDURE  DEFINITION  OF  SetLo 

This  sets  the  global  variable  Lo  after  y  and  A  have  been  set  (see  Algorithm  II. 
local  constants  required  for  SetLo: 
k  s  wa-Bo 

implementation  of  SetLo: 

U  =  [k/lK-A)]’'"'*”  Equation  7 

PROCEDURE  DEFINITION  OF  SetLo&A(L,  s) 

This  sets  the  gobal  variables  A  and  Lg  when  y  is  known  (see  Algorithm  1)  so 
that  ^0  includes  the  point  indicated  by  L  and  s.  A  binary  search  is  used. 
local  constants  required  for  SetLo&A(L,  si: 
k  s  u^-a-Bo 

local  variables  required  for  SetLo&A(L,  s): 

AA,  AL 

implementation  of  SetLo&A(L,  s): 

A  =  k/(K'L’’*M  From  Equation  7;  the  maximum  value  of  A  is  found 

SetLo  when  L  is  assumed  to  be  Lg  at  dusk. 

AL  =  L  -  L^(s)  from  Algorithm  4. 

A  A  =  A  This  is  the  starting  step  size  of  search. 

while  abs(AL)  >  ~  repeat  the  following  loop: 
begin  while  loop 

A  A  =  Vi  A  A  The  step  size  is  halved  during  each  iteration. 

if  AL  <  0  then 

A  =  A  +  AA  A  was  too  small. 

else 

A  =  A-AA  A  was  too  large. 

SetLo 

AL  =  L  -  L^ls)  from  Algorithm  4. 

end  while  loop 

Eq.  (3).  For  a  given  L-shell  and  local  time,  y  and  W  are  first  determined  from  Eqs.  (3) 
and  (4)  using  the  above  process,  and  then  U  is  determined  using  Eq.  (5).  The  given  L- 
shell  is  a  valid  solution  if  the  difference,  A(7,  between  Uy  and  U  is  zero.  Algorithm  6 
uses  this  process  to  solve  for  AL/  for  a  measured  position  and  energy.  Valid  solutions 
of  L  for  each  local  time  give  the  drift  path  of  the  particle. 

In  summary,  the  drift  paths  of  a  particle  are  determined  by  calculating  the  three 
constants  of  motion,  /i„,  for  its  kinetic  energy  and  position  in  space  using  Eqs. 

(3),  (4),  and  (S).  Using  these  constants  of  motion,  Eqs.  (3)  and  (4)  are  solved 
simultaneously  to  give  kinetic  energy  as  a  function  of  L-shell.  All  L-shells  at  a  given 
local  time  for  which  the  difference  between  Un  and  U  [Eq.  (5)]  is  zero  are  on  the  drift 
path  of  the  particle. 


4,  CALCULATION  OF  ELECTRON  CUTOFF  ENERGY,  Wc 

In  a  manner  similar  to  that  used  to  calculate  drift  paths,  the  range  in  energy  of 
electrons  coming  from  the  tail  that  should  be  observed  at  a  given  position  in  the 
magnetosphere  can  be  determined  using  Eqs.  (3),  (4),  and  (S).  It  is  known  that  electrons 
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Figure!  Exan^Ieofthe  iterative  solution  ofEquation  9  for  1.69.  Solution  starts  with  Point 
A  and  converges  to  Point  D.  Each  thin  line  represents  a  linear  approximation  to  I(y)  at  one  of 
the  points. 


with  W>  0  execute  a  gradient  curvature  drift  in  the  direction  of  the  corotation  drift.  For 
increasing  kinetic  energy  this  gradient  curvature  drift  comes  to  dominate  the  electric  field 
drift,  and  since  this  drift  has  no  radial  component,  the  electron  trajectory  becomes  more 
circular  and  eventually  closes.  For  any  position  outside  of  electrons  measured  with 
W=0  are  on  open  drift  paths  as  shown  in  Figure  1,  and  electrons  measured  with  W-+oo 
are  on  closed  circular  drift  paths  because  the  magnetic  drift  completely  dominates  the 
electric  field  drift.  For  the  position  where  the  electron  is  observed,  ftere  is  a  cutoff 
energy,  Wc,  at  which  the  drift  trajectory  changes  from  open  to  closed.  Electrons 
convecting  from  the  tail  with  energies  ranging  from  0  to  Wc  should  be  observed  at  that 
position.  Wf.  is  not  defined  for  positions  inside  where  all  electrons  move  on  closed 
drift  paths. 

Figure  3  illustrates  the  electron  cutoff  energy,  W^.  Six  sets  of  electron  drift  paths 
in  the  equatorial  plane  are  shown  in  this  figure,  and  these  are  calculated  using  the  method 
present^  in  Section  3  and  Algorithm  6.  The  constants  of  motion  are  determined  for 
electrons  with  one  of  six  different  values  of  measured  at  the  white  circle  with 
Lj^=6.0  Re,  yji/=0.90,  and  sin(^j^)=-0.95.  The  six  values  of  are  given  on  the  right 
side  of  the  figure,  and  are  labeled  A  to  F.  In  this  example  Lo-5  Re,  and  the  shaded 
region  is  bound  by  Xq.  Line  A  shows  the  path  for  electrons  with  Wj^=0.  In  this  case 
there  is  a  closed  eastward  path  around  the  Earth  inside  of  Xo,  and  there  is  a  sq)arate 
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Algorithm  6  Calculation  of  AU  for  given  constants  of  motion  and  position. 


FUNCTION  DEFINITION  OF  AiUfj.  J.  U,  electron.  L,  s) 
local  constants  required  by  AUOi/,  J,  U,  electron,  L,  s): 

k  Si  Bo-a*-w 

local  variables  required  for  AU(//,  J,  U,  electron,  L,  s): 

W,  y  kinetic  energy  (eV)  and  sine  of  eguatoriai  pitch  angle 

a,,  b,  linear  coefficients  of  Ij 

i,  j  integer  indices  of  4 

Q  of  Equation  9. 

implementation  of  AU(//,  J,  U,  electron,  L,  s): 
if  p  =  0.0  then 

W  =  <J/L/lo)*  Equation  4. 

else  if  J  =  0.0  then 

W  =  pIL^  Equation  3. 

else 


begin  conditional  block 

i  =  N  N  is  from  Algorithm  1. 

j  =  0  this  ensures  that  i  /  the  first  time. 

Q  =  3- (Up)*  Equation  9. 

while  i  j  repeat  the  following  loop:  Check  for  convergence 


begin  while  loop 

j  =  i 

b.  =  N-|l,-l,,) 
a,  =  I.  -  b,-i/N 
y  =  a,  /  IQ-b,) 
i  =  trunclN-y)  +  1 
end  while  loop 
W  =  M<yL)*U 
end  conditional  block 
if  electron  then 

return  value  of  U  -  W  +  a- 

else 


slope  between  4  and  4.,. 
zero  offset. 

Equation  10. 

update  index  with  new  value  of  y. 
Equation  3. 


•s-L*'  -  k/L  Equation  5,  q=^-1 


return  value  of  U  -  W  -  a-A-s-L*'  +  k/L  Equation  5,  q  =  +  1 


sunward  path  on  the  dusk  side  of  the  Earth  that  is  open  to  the  tail.  The  point  of 
measurement  is  on  the  open  path.  Line  B  shows  the  path  for  H^=0.6  keV  electrons. 
Here  the  inner  and  outer  paths  join  at  a  single  stagnation  point  on  the  dusk  meridian  near 
5.2  Re,  and  the  point  of  measurement  is  located  on  the  open  outer  open  path.  Line  C 
shows  paths  for  Wj^=3.8  keV  electrons.  Because  the  inner  and  outer  paths  have  joined, 
only  one  path  exists,  and  this  is  open  to  the  tail  and  crosses  the  dawn-dusk  meridian  on 
the  dawn  side  of  the  Earth.  Line  D  shows  the  path  for  electrons  with  y^u—^c-  Here 
a  closed  eastward  path  joins  at  a  single  point  on  the  dusk  meridian  near  6.5  Re  with  an 
open  sunward  path.  This  path  looks  very  similar  to  that  of  Line  B,  but  here  the  point  or 
measurement  is  on  the  closed  inner  path  of  Line  D.  Line  D  shows  the  transition 
trajectory  between  open  and  closed  drift  paths,  because  all  electrons  with  W^<  are 
on  open  paths  while  higher  energies,  as  shown  by  Lines  E  and  F,  place  the  point  of 
measurement  on  closed  eastward  paths.  In  the  last  two  cases,  the  open  paths  from  the 
tail  and  the  closed  eastward  paths  are  completely  separate. 


To  calculate  Wc  for  a  measured  location,  the  value  of  W],,  is  determined  that 
yields  a  single  solution  of  AU-0  on  the  dusk  meridian  as  shown  by  Line  D  in  Figure  3. 


10 


shaded  region  is  bounded  by  and  arrows  show  the  direction  of  drift  along  the  paths. 

In  Figure  4,  Algorithm  6  is  used  to  calculate  6XJ,  and  each  curve  shows  MJ  versus  L- 
shell  at  dusk  for  the  same  six  values  of  and  for  the  same  measured  location  given  for 
Figure  3.  The  peak  values  of  At7  at  dusk  are  referred  to  as  £^Vp.  For  Curves  A,  E,  aiMl 
F,  AUp>0,  so  there  are  two  solutions  to  AU=0  at  dusk  for  these  energies.  Because 
there  are  two  solutions,  two  separate  paths  are  seen  at  dusk  in  Figure  3  for  Lines  A,  E, 
and  F.  For  Curve  C,  AUp<  0,  so  that  there  are  no  solutions  to  A17=0,  and  no  paths  are 
seen  at  dusk  in  Figure  3  for  Line  C.  For  Curves  B  and  D,  AUp=0,  so  that  both  have 
exactly  one  solution  for  AU=0  at  dusk.  Both  Lines  B  and  D  in  Figure  3  have  exactly 
one  path  on  the  dusk  meridian,  and  their  respective  inner  and  outer  paths  join  at  this 
point.  Curve  D  shows  the  proper  solution  for  W^.,  as  was  discussed  in  the  previous 
paragraph,  because  it  places  the  point  of  measurement  on  the  inner  path,  which  is  closed 
for  W^>  Wc  and  open  for  W^<  When  two  values  of  exist  for  which  AUp=0, 
W(~  is  equal  to  the  greater  value. 
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4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0  B.5  9.0 


L-shell  at  dusk.  Re 

F^ure  4  Al/  as  a  function  of  L-shell  at  dusk  (sin(^)=-l)  for  the  same  conditions,  measured 
location,  and  energies  as  those  given  for  Figure  3. 

In  order  to  determine  the  value  of  the  L-shell  at  dusk  at  which  the  peak 
occurs  is  identified,  and  this  L-shell  is  referred  to  as  Lp,  so  that  AUp=AU(Lp).  occurs 
at  the  L-shell  where  dAUIdL=0  at  dusk.  AIJ=^U^U,  and  U  comes  from  Eq.  (5)  with 
9=-l  and  sin(<^)=-l,  so  that 

=  0  (11) 

oL  oL 

The  term  dWIdL  is  derived  from  the  simultaneous  solution  of  Eqs.  (3)  and  (4)  as  is  done 
in  Algorithm  6.  From  Eq.  (3), 

^  (12) 

The  dyIdL  term  can  be  derived  from  Eq.  (10)  as 

I  -  (*3) 


and  the  dQ/dL  term  can  be  derived  from  Eq.  (9)  as 


Algorithm  7  Calculation  of  SW/dL  as  a  function  of  L'Sheli  at  dusk  for  a  given  set  of  constants 
of  motion. 


FUNCTION  DEFINITION  OF  aW/aL(p.  J,  L) 
local  variables  required  for  dWIBUjj,  J,  L): 

Dw  value  of  dW/dL  returned  by  function. 

j  integer  indices  to  If. 

a„  b,  linear  coefficients  to  If 

Q,  Dq,  y,  Dy 

implementation  of  SWIdLip,  J.  L): 
if  //  =  0  then 

Dw  =  -2[J/(L*lo)l^/L  derived  from  Equation  4. 
else  if  J  =  0  then 

Dw  =  -2plL*  derived  from  Equation  3. 

else 


begin  conditional  block 
i  =  N  N  is  from  Algorithm  1. 

j  =  0  this  ensures  that  i-Aj  the  first  time. 

Q  =  J’lL///)**  Equation  9. 

Dq  =  V4Q/L  Equation  14. 

while  i  ^  i  repeat  the  following  loop: 
begin  while  loop 


j  =  i 

b,  =  N  (l,  - 1,,) 
a,  =  1;  -  b,-i/N 
V  =  a,  /  (0  -  b.) 
i  =  truncly-N)  +  1 
end  while  loop 
D,  =  -a,-DQ  /  (Q  -  b,)* 

Dw  =  -;/-y^-L-*(2D/L  +  3 
end  conditional  block 
return  the  value  of  Dw 


slope  between  4  and 
zero  offset. 

Equation  10. 

update  index  with  new  value  of  y. 

Equation  13. 
y)  Equation  12. 


4?  -  'ACL-*  (14) 

oL 

TTie  solution  to  dWIBL  is  implemented  in  Algorithm  7.  If  either  ^  or  /  is  equal  to  zero 
then  dWIdL  is  solved  from  either  Eq.  (4)  or  Eq.  (3).  Otherwise,  the  iterative  search 
shown  for  Algorithm  6  is  used  to  solve  for  y  and  W,  and  Eqs.  (12),  (13),  and  (14)  are 
used  to  solve  for  dWIdL. 

There  is  no  analytic  solution  to  Eq.  (1 1)  when  Eq.  (14)  is  substituted  for  dWIdL, 
so  a  numeric  solution  is  used.  First  dWIdL  is  approximated  by  a  simple  form  which 
gives  an  analytic  solution  to  Eq.  (11).  Inspection  of  the  two  extreme  cases,  Pu—O  and 
Jm—0,  of  Algorithm  6  shows  that  dWIdL  «  kL^  where  m  ranges  from  3  to  4.  This  is 
approximated  by  the  form 

«  +  (IS) 

quadratic  equation  when  it  is  substituted  into  Eq. 


dW 

"dL 

This  form  is  used  because  it  yields  a 
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Algorithm  8  Calculation  of  Lp  for  givoi  values  of  n  and  J  for  electrons. 


FUNCTION  DEFINITION  OF  PeakWcl//.  J) 
local  constants  required  for  PeakWcO^.  J): 

k  3s  Bo-a*-a; 

local  variables  required  for  PeakW^ijj,  J): 

L,  L,,  Lj  These  hold  the  current  solutions  of  L-shell 

D,  D,,  Dj  These  hold  dW/dL  for  the  current  solutions  of  L-shell 

c,,  C2»  Qa/  Qb*  Qq 
implementation  of  PeakWclp,  J): 

L,  =  Lo 

D,  =  3W/3LU/.  J,  L,) 

Lj  =  L,  +  1.0 
Dj  =  3W/3L(;/,  J,  Lj) 

Qa  =  -ra-A 

while  absIL^'L,)  >  3,.  repeat  the  following  loop: 
begin  while  loop 
Cj  =  (L,2-D,  - 

c,  =  L,*‘Di  -  Cj'L, 

Qg  =  k  -  c^ 


Find  two  starting  points  for  the  iteration. 

Lg  is  usually  within  afew  of  the  final  solution. 


First  coefficient  of  Equation  -1 6. 


Qq  “  Qe  +  4'Q^'C2 
L  =  [ysl-Qg  -  Qq’‘)/Qa1’'"'"” 

D  =  awidup,  J.  Li 
if  abs(L  *  L,)  <  abs(L  >  Lj)  tfien 
begin  conditional  block 
L2  =  L 
Dj  =  D 

end  conditional  block 

else 

begin  conditional  block 
L,  =  L 
D,  =  D 

end  conditional  block 
end  while  loop 
return  value  of  L 


Solution  to  Equation  15 
coefficients. 

Second  coefficient  of  Equation  1 6. 
Value  inside  of  root  in  Equation  1 7. 
Equation  1 7  yields  better  value  for  L. 
New  value  for  5W/bL  is  computed. 

Replace  whichever  solution  is 
farthest  from  new  solution. 


(11).  Substitution  of  Eq.  (15)  into  Eq.  (1 1),  and  multiplying  both  sides  by  gives 

0  »  -  (fiyw-c,)Ly^  -  c, 


The  solution  of  interest  is 


22. 


(17) 


where  QA~-yaA,  Qb-B^u-c^,  and  Qc^-c^.  The  values  of  and  Cj  are  first  determined 
from  two  points  at  L=Lg  and  L=Lg-\-\.  An  estimate  of  Lp  is  determined  using  Eq.  (17), 
and  this  new  value  is  used  to  replace  the  initial  estimated  value  of  L-shell  farthest  from 
the  new  estimate.  New  values  for  q  and  Cj  are  determined  which  then  give  a  better 
estimate  of  Lp.  This  process  is  repeated  until  it  converges  to  a  solution,  usually  in  four 
to  five  iterations.  Algorithm  8  shows  this  method  of  determining  Lp. 
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L-shell  at  dusk.  Re 


Figure  5  Slope  of  Curve  D  in  Figure  4.  The  bold  line  is  the  actual  slope,  and  the  thin  solid  lines 
are  approximations  using  Eq.  (IS).  The  dashed  lines  exaggerate  the  error  betwe«i  the 
approximations  and  actual  slope  by  a  factor  of  SO. 

Figure  5  illustrates  Algorithm  8  for  the  solution  of  Lp  for  the  case  W^, 
Curve  D  of  Figure  4.  The  bold  line  shows  d^UldL  as  a  function  of  L-shell  at  dusk.  Lp 
is  located  where  this  line  crosses  zero  (Point  T).  Points  P  and  Q  at  L=Lo  and  L=Z,+ 1 
are  used  to  determine  initial  values  for  c,  and  [Eq.  (15)].  Line  u,  which  is  a  thin  solid 
line,  shows  the  approximation  to  dAUIdL  by  using  Eq.  (15).  Because  the  approximation 
is  very  close  to  the  actual  solution,  Line  u',  which  is  a  thin  dashed  line,  shows  the  error 
between  the  approximation  and  the  actual  solution  exaggerated  by  a  factor  of  50. 
Equation  (17)  identifies  the  L-shell  of  the  zero  intersect  of  Line  u  at  Point  R.  The  second 
iteration  uses  Points  Q  and  R  to  determine  c,  and  Cj.  This  approximation  is  shown  by 
Line  v  with  the  exaggerated  error  shown  by  Line  v'.  Here  Ae  approximation  is  more 
accurate  near  the  solution.  Point  S,  which  is  the  zero  intercept  of  Line  v,  is  used  with 
Point  R  for  the  final  iteration  shown  by  Line  w.  Line  w'  shows  that  the  approximation 
is  very  good  near  the  solution  shown  by  Point  T.  The  iteration  is  stopped  here  because 
the  difference  between  the  Points  S  and  T  is  sufficiently  small. 

Because  there  are  potentially  two  values  of  W„  for  which  AUp=Q,  as  shown  by 
Curves  B  and  D  in  Figure  4,  it  is  necessary  to  identify  a  lower  limit  to  which  is 
always  larger  than  the  smaller  solution  to  AUp=0  and  is  less  than  or  equal  to  the  larger, 
correct  solution.  In  Figure  4,  all  of  the  curves  have  the  same  value  of  AU  at  L=6.0  Re 
because  this  is  the  value  of  L^.  When  L^L^,  conservation  of  the  first  two  adiabatic 
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invariants  [Eqs.  (3)  and  (4)]  requires  that  and  y=yu.  Substituting  Eq.  (5)  into 

AU=Uj^-U  yields 

AU(L=LJ  =  -aALl[\  +sin(<^J]  (18) 

and  this  shows  that  AI7(L=Lj^)  <0.  As  can  be  seen  by  looking  at  Curve  C,  AU(L=Lm) 
is  the  minimum  possible  value  of  AUp.  The  value  of  that  gives  AUp-AU{L=Ly)  is 
referred  to  as  and  is  obtained  by  setting  dAUIdLiL=L^)  to  zero.  Since  is  a 
constant,  the  equation  is 

^(L=LJ  =  -yaALl'^  =  0  (19) 

For  two  different  electrons,  labeled  1  and  2,  if  and  y,=y2  then  Eqs.  (3)  and  (4) 
require  that respectively.  With  these  requirements,  Eq.  (9) 
shows  Q2~Git  Eq.  (14)  shows  hQ^dL=dQjdLy  and  Eq.  (13)  shows  dy2ldL-dy^ldL. 
Equation  (12)  then  shows  that 


aw, 

W,  =  W, — I 
'  ‘  dL 


aw, 

~dL 


(20) 


To  solve  for  W^  dWldL(L=Lu)  from  Eq.  (19)  is  substituted  into  dW^ldL  in  Eq.  (20). 
Any  energy,  W^,  can  be  used  for  W,  provided  that  dW^ldL  uses  the  constants  of  motion 
from  L=Lu,  y=yu,  <l>=4>u,  and  W=  W^.  This  gives 

Substitution  of  Eq.  (7)  into  Eq.  (21)  shows  that  W<,>0  for  L„>L^,  and  Wo<0  for 
Lu<Lq.  When  Wo>0,  two  values  of  may  exist  for  which  A(/,=0,  and  W„  is 
between  these  two  values.  In  Figure  4,  W,,  is  the  energy  used  for  Curve  C,  and  Curves 
B  and  D  use  the  lower  and  upper  values  of  W^  for  which  AUp=Q.  When  W^^O,  only 
one  solution  can  exist  because  kinetic  energy  cannot  be  negative. 

Wc  is  determined  by  identifying  the  value  of  W^>Wp  that  gives  AUp=0. 
Algorithm  9  shows  the  implementation  of  the  iterative  search  used  to  find  W^.  for  a 
position  given  as  y^yu,  L-L^,  and  sin(<^)=S;^.  If  Ly<L^,  then  the  location  L^f,  is 
outside  of  Xo  if  U  for  Wj,,=0  electrons  is  less  than  U  of  This  can  be  seen  in  Figure 
1  in  which  the  contour  lines  show  negative  U  because  kinetic  energy  is  zero  and  q  is 
negative.  If  the  measured  position  is  inside  of  Xo,  then  W^.  is  given  the  value  of  -1  to 
indicate  that  is  not  defined.  The  initial  lower  value  of  W„  for  the  search  is  given  by 
Eq.  (21).  This  equation  gives  a  negative  solution  if  and  for  these  cases  an  initii 

lower  value  of  W„=0  is  used.  The  initial  upper  value  is  chosen  arbitrarily  to  be  10  keV 
greater  than  the  lower  value.  A  third  point  is  picked  by  assuming  the  first  two  points  lie 
on  a  parabola,  and  the  slope  of  the  parabola  is  zero  at  the  lower  point.  If  W0<  0,  a  third 
point  is  simply  picked  halfway  between  the  first  two  points.  A  parabola  is  then  used  to 


dW^ 

1l' 


(W^WX^L^) 


“1 


(21) 
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Algorithm  9  Calculation  of  electron  cutoff  energy,  W^,  for  a  position  given  by  L,  y,  and  sin(^). 


FUNCTION  DEFINITION  OF  Wc(y.  L,  s) 
local  constants  required  for  Wcly,  L,  s}: 

k  5s  Bo'wa* 

local  variables  required  for  Wc(y.  L,  s): 

Uo.  n,.  Wo,  W,.  Wj,  AUo,  AU„  AUj.  Q*,  0^.  Qc 
implementation  of  Wgly,  L,  s): 

=  p(y.  1 ,  U  Value  of  at  y,L,s  if  W=  1  eV;  iJ  =  fi,-W  Eqn  3 

J,  =  J{y,  1 .  U  Value  ofJat  /,L,s  if  W=  1  eV;  J  ^  J,-W^  Eqn  4 

Uo  =  U(true,  0,  L,  s)  Value  of  U  at  y,L,s  if  W=0  eV;  U  =  Ug  +  W  Eqn  5 
if  (L  <  Lo)  and  [AU(0,  0,  Uo,  true,  Lo.  -1)  >  01  then  satellite  is  inside  of  Xg 
return  value  of  -1  and  exit  function  -1  indicates  that  Wg  does  not  exist 
Equation  2 1  finds  initial  starting  point  of  search. 
if  L  >  Lo  then 

Wo  =  -  ya-A-L*^’)  /aW/3L|;/,  W„^,  J,  W^*,  L) 

else 

Wo  =  0  Equation  21  returns  a  negative  value  so  use  Wg  =  0 
AUo  =  AUU/rWo,  J.  Wo**,  Uo+Wo,  true,  PeakWe07,  Wo,  J.  Wo'^),  -1)  * 

The  second  starting  point  is  offset  from  the  first  by 
W,  =  Wo  +  W„^ 

AU,  =  AU(/i,  W„  Ji-W,**,  Uo  +  W„  true,  PeakWcU/,  W„  -1)  * 

A  parabolic  approximation  is  used  to  find  the  third  point.  It  uses  the  points  Wg 
and  W„  and  assumes  the  slope  at  Wg  is  zero. 

If  Wo  =  0  then 

Wj  =  %  W,  Slope  at  Wg  is  not  zero  if  Wg=0;  use  midpoint. 

else 

Wj  =  Wg  +  W^TAUo/IAUo-AU,)]*  This  is  zero  intercept  of  parabola. 

AU2  =  AU(//,  W2,  Ji-W/,  Uo  +  Wj,  true,  PeakWcU/rWj,  J.-Wz**),  -1)  * 
Reorder  W„  and  Wj  so  that  |A6/o|  <  |Ay,|  <  |  Af/j| 
if  abs(AU,)  <  abslAUol  then 

swap  values  of  AU,  with  AUo  s^d  W,  with  Wo 
if  abslAUj)  <  abslAU,)  then 

swap  values  of  AU2  with  AU,  and  W;  with  W, 
if  abslAU,)  <  abs(AUo)  then 

swap  values  of  AU,  with  AUo  and  W,  with  Wg 
while  abs(Wi-Wo)  >  dw'<Wo  +  W,)  repeat  following  loop: 
begin  while  loop 

Q*  =  lAUo  lWj-W,)  +  AU,  (Wo-W2)  +  AU2  (W,-Wo)]  / 

[Wo^-lWj-W,)  +  W,2-<Wo-W2)  +  W2*-(W,-Wo)] 

Qe  =  lAUo-AU,  -  0*MWo"-W,*)]  /  (Wg-W,) 

Qc  =  AUo  -  Wo  (Qb  +  QvWg) 

AUj  *  AU,  Shift  point  1  to  point  2 

Wj  =  W,  and  point  0  to  point  1, 

AU,  =  AUo  and  find  new  point  0. 

W,  =  Wg 

Wo  =  ’A-[-Qe  +  (Qfl*  -  A-Q^Qc)*]  /  Qa  Quadratic  Formula 

AUo  =  AUU/,-Wo,  J,-Wo*,  Uo  +  Wg,  true,  PeakWc(it/,  Wo,  J,  Wo'‘),-1 )  * 
end  while  loop 
return  value  of  Wo 

★  These  lines  can  be  made  more  efficient  by  calculating  p,-W„  and  J,-W„’^  only  once 
for  each  line  instead  of  twice. 
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Figure  6  Relationship  between  AUp  and  for  the  conditions  given  for  Figure  4.  This 
illustrates  the  solution  for 


approximate  the  relationship  between  AUp  and  using  the  three  points,  and  it  is  used 
to  solve  for  a  new  value  of  The  point  that  is  farthest  from  AUp=0  is  replaced  by 
the  new  solution  for  in  the  next  iteration.  The  iteration  is  continued  until  the  value 
of  is  known  to  sufficient  accuracy,  usually  in  about  three  iterations. 

Figure  6  illustrates  the  solution  for  Wf.  for  the  conditions  given  for  Figures  3  and 
4,  and  Points  A  through  F  have  the  same  energy  as  do  those  in  the  previous  figures.  The 
X  axis  shows  the  y  axis  shows  AUp,  and  the  dark  curve  shows  the  relationship 
between  the  two.  The  initial  lower  boundary  at  Point  C  is  identified  using  Eq.  (21),  and 
die  initial  upper  boundary  at  Point  F  is  found  by  adding  10  keV  to  the  energy  of  Point 
C.  Point  E  is  found  using  the  zero  intercept  of  the  parabola  formed  from  Points  C  and 
F  with  a  slope  of  zero  at  Point  C  (shown  by  thin  dotted  line).  The  zero  intercept  of  the 
parabola  formed  from  Points  C,  E  and  F  (thin  solid  line)  is  very  close  to  the  true  solution 
shown  by  Point  D.  The  process  converges  to  Point  D  on  the  next  iteiation. 


5.  CALCULATION  OF  ION  LOWER  CUTOFF  ENERGY,  Wi 

Ions  move  in  a  more  complex  manner  than  do  electrons  because  the  magnetic 
gradient  and  curvature  drift  of  ions  is  westward  which  is  counter  ta  that  of  the  corotation 
electric  field  drift.  Ions  with  zero  kinetic  energy  experience  no  gradient  and  curvature 
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drift,  and  so  move  identically  to  electrons  with  zero  kinetic  energy.  Therefore,  an  ion 
inside  of  with  W=0  moves  eastward  on  a  closed  drift  path.  As  W-*oo,  the  gradient 
curvature  drift  completely  dominates  the  electric  field  drift  so  that  the  ion  moves 
westward  on  a  closed  circular  drift  path.  The  drift  paths  for  intermediate  kinetic  energies 
transition  ft-om  closed  eastward  paths  to  closed  westward  paths,  and  many  of  these  paths 
are  open  to  the  tail.  The  bounding  limits  to  the  range  of  kinetic  energy  open  to  the  tail 
are  for  the  lower  bound  and  W„  for  the  upper  bound.  Ions  with  W<  move  on 
closed  eastward  drift  paths,  and  ions  with  W>  move  on  closed  westward  drift  paths. 

is  not  defined  for  positions  outside  of  because  here  no  ions  can  move  on  closed 
eastward  drift  paths.  The  two  ion  cutoff  energies  are  discussed  separately,  in  this 
section,  and  in  Section  6. 

Figure  7  illustrates  the  ion  lower  cutoff  energy,  Four  sets  of  ion  drift  paths 
in  the  equatorial  plane  are  shown  in  this  figure;  they  are  calculated  using  the  method 


K,  O.OkeV 
B,  1.0keV=W, 


X  ^  ,  - - - -  7.2keV 

00  D,  9.0keV 

Figure  7  Convection  paths  of  singly  charged  ions  measured  at  the  white  circle  for  four  values 
of  Wu.  The  shaded  region  is  bounded  by  ifo.  arrows  show  the  direction  of  drift  along  the 


presented  in  Section  3  and  Algorithm  6.  The  constants  of  motion  are  determined  for  ions 
with  one  of  four  different  values  of  measured  at  the  white  circle  with  Lj^—3.5  Rg, 
>>^=0.90,  and  sin(<^v)=0-95.  In  each  case,  W^<  W„.  The  four  values  of  are  given 
on  the  right  side  of  the  figure  and  are  labeled  A  to  D.  In  this  example  Lq=6  Rg,  and  the 
shaded  region  is  bounded  by  iCo-  Line  A  shows  the  path  for  ions  with  Wj|^=0.  For  these 
ions,  there  is  a  closed  eastward  path  around  the  Earth  inside  of  Xq,  and  there  is  a 
separate  sunward  path  that  is  open  to  the  tail  on  the  dusk  side  of  the  Earth.  The  point 
of  measurement  is  on  the  clos^  path  whereas  it  is  on  the  open  path  for  electrons  with 
WJ^,=0.  Line  B  shows  the  path  for  ions  with  Wj^=  Here  a  closed  eastward  path  joins 
at  a  single  point  on  the  dusk  meridian  with  an  open  sunward  path.  This  path  looks  very 
similar  to  that  for  electrons  with  except  that  the  paths  Join  inside  of  Xq.  This 

is  the  transition  trajectory  between  open  and  closed  drift  paths,  because  all  ions  with 
are  on  closed  paths  while  higher  energies,  as  shown  by  Lines  C  and  D,  place 
the  point  of  measurement  on  paths  open  to  the  tail.  In  the  last  two  cases,  the  open  paths 
from  the  tail  and  the  closed  eastward  paths  join  to  form  a  single  path  that  crosses  the 
dawn-dusk  meridian  on  the  dawn  side  of  the  Earth.  Lines  C  and  D  also  show  two  closed 
westward  drift  paths  inside  of  the  point  of  measurement.  Electron  drift  paths  have  no 
equivalent  to  these  paths  because  the  electron  gradient  curvature  drift  is  eastward. 

is  calculated  in  a  manner  similar  to  that  of  W^.  In  Figure  8,  Algorithm  6  is 
used  to  calculate  AU,  and  four  curves  show  AU  versus  L-shell  at  dusk  for  the  same  four 
values  of  and  for  the  same  measured  location  given  for  Figure  7.  Comparison  of  this 


Figure  8  AU  as  a  function  of  L-shell  at  dusk  (sin(^)=-I)  for  the  same  conditions,  measured 
location,  and  energies  as  those  given  for  Figure  7.  The  crosses  mark  the  location  of  AUp. 
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figure  to  Figure  4  shows  the  more  complex  nature  of  the  ions.  In  Figure  4  there  is  one 
peak  for  each  electron  energy.  In  Figure  8,  peaks  can  be  seen  in  Curves  A,  B,  and  C 
at  the  crosses,  but  the  peaks  are  negative,  and  a  peak  does  not  exist  for  Curve  D.  Also, 
Curves  C  and  D  show  an  extra  solution  to  AU—0  at  L-shells  between  2  and  3  Rg.  These 
two  inner  solutions  form  the  closed  westward  drift  paths  for  Lines  C  and  D  of  Figure  7. 
The  peak  values  of  AU  at  dusk  are  referred  to  as  AUp.  In  cases  where  a  peak  does  not 
exist,  AUp  is  defined  at  the  point  were  the  slope  minimizes  as  is  shown  by  the  cross  on 
Curve  D.  For  Curve  A,  AUp<0,  so  there  are  two  solutions  to  AC/=0  at  dusk  for  this 
energy.  Because  there  are  two  solutions,  there  are  two  separate  paths  seen  at  dusk  in 
Figure  7  for  Line  A.  For  Curves  C  and  D,  AUp>0,  so  that  there  are  no  solutions  other 
than  that  for  the  inner  closed  westward  drift  path,  and  only  this  path  is  seen  at  dusk  in 
Figure  7  for  Lines  C  and  D.  AUp—Q  for  Curve  B,  so  Line  B  of  Figure  7  has  exactly  one 


Algorithm  10  Calculation  of  Lp  at  dusk  for  given  values  of  and  J  for  ions. 

FUNCTION  DEFINITION  OF  PeakW^O/,  J) 
local  constants  required  for  PeakW^I/i,  J); 
k  =  -Bo-a^n; 

local  variables  required  for  PeakV/^ifj,  J): 

L,  L,,  Lj  These  hold  the  current  solutions  of  L-shell 

D,  D,,  Dj  These  hold  dW/dL  for  the  current  solutions  of  L-shell 

c,,  C2,  Q-p,  Qg,  Qq 
implementation  of  PeakW,^!;/,  J): 

L,  =  Lq  Find  two  starting  points  for  the  iteration. 

D,  =  3W/3L0:/,  J,  L,)  Lgis  usually  within  a  few  Rgof  the  final  solution. 

L2  =  L,  -  1 .0 

Dj  =  3W/3L(;/,  J,  Lj) 

=  K  a-A  First  coefficient  for  Equation  23. 

while  abslLj-L,)  >  repeat  the  following  loop: 
begin  white  loop 

C2  =  IL,^-D,  -  L2^'D2)  /  (L/''’  -  L2'"’>  Solution  to  Equation  15 
c,  =  Li^  D,  -  C2-L,  ’^’  coefficients. 

Qg  =  k  -  c,  Second  coefficient  for  Equation  23. 

Qq  =  Qb^  +  4  Qa-C2  Value  inside  of  root  in  Equation  23. 

if  Qq  <  0  then  If  no  peak  is  found  then  search  for 

Qq  =  0  minimum  slope. 

L  =  [’AI-Qb  +  Qq’‘)/QJ"'’'*”  Equation  23  yields  better  value  for  L. 

D  =  3W/3L{//,  J,  L)  New  value  for  dW/dL  is  computed. 

if  absIL  -  L,)  <  abs(L  -  L2)  then  Replace  whichever  solution  is 

begin  conditional  block  farthest  from  new  solution. 

L2  =  L 
□2  =  0 

end  conditional  block 

else 

begin  conditional  block 
L,  =  L 
D,  =  D 

end  conditional  block 
end  while  loop 
return  value  of  L 
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solution  for  AC/=0  at  dusk,  and  the  inner  and  outer  paths  join  at  this  point. 


To  determine  the  value  of  At/„  the  L-shell  at  dusk  at  which  the  peak  occurs  is 
identified.  TTiis  L-shell  is  referred  to  as  Lp,  so  that  AUp-AU(Lp).  Lp  occurs  at  the  L- 
shell  where  dAUIdL—0  at  dusk.  AU=  Uu-U,  and  U  comes  from  Eq.  (5)  with  q—  1  and 
sin(<#>)=-l,  so  that 

^(L=L,)  =  -^(L=Lp)^yaALr'  =  0  (22) 

uL  OL 

This  solution  is  identical  to  Eq.  (11)  except  for  sign  changes  in  the  second  and  third 
terms.  It  is  solved  iteratively  as  is  Eq.  (11)  by  using  Eq.  (15)  except  that  the  solution 
of  interest  is 


-QsHQl-^Q^Qcr 

2G/ 


(23) 


where  Q^=yaA,  QB=-Bf/a^u-Ci,  and  Qc=-C2-  The  values  of  c,  and  Cj  of  Eq.  (15)  are 
first  determined  from  initial  estimates  of  L=Lo  and  L=Lo-l.  An  estimate  of  Lp  is 
determined  using  Eq.  (23),  and  this  new  value  is  used  to  replace  the  initial  estimated 
value  of  L-shell  farthest  from  the  new  estimate.  New  values  for  Ci  and  Cj  are  determined 
giving  a  better  estimate  of  Lp.  If  no  peak  is  found,  such  as  in  the  case  shown  by  Curve 
D  of  Figure  8,  the  term  in  Eq.  (23)  inside  of  the  radical  is  set  to  zero,  so  that  Eq.  (23) 
identifies  the  L-shell  of  minimum  slope.  This  process  is  repeated  until  it  converges  to 
a  solution,  usually  in  four  to  five  iterations.  Algorithm  10  shows  this  method  of 
determining  Lp. 

Figure  9  illustrates  Algorithm  10  for  the  solution  of  Lp  for  the  case 
Curve  B  of  Figure  8.  The  bold  line  shows  dAUIdL  as  a  function  of  L-shell  at  dusk.  Lp 
is  located  where  this  line  crosses  zero  (Point  S).  Points  P  and  Q  at  £.=Lo  and  L=Io-l 
are  used  to  determine  initial  values  for  c,  and  Cj  [Eq.  (15)].  Line  u,  which  is  a  thin  solid 
line,  shows  the  approximation  to  dAUIdL  using  Eq.  (15).  Because  the  approximation  is 
very  close  to  the  actual  solution.  Line  u',  which  is  a  thin  dashed  line,  is  shown  with  the 
error  between  the  approximation  and  the  actual  solution  exaggerated  by  a  factor  of  500. 
Equation  (23)  identifies  the  L-shell  of  the  zero  intersect  of  Line  u  at  Point  R.  The  second 
and  final  iteration  uses  Points  P  and  R  to  determine  q  and  Cj.  This  approximation  is 
shown  by  Line  v  with  the  exaggerated  error  shown  by  Line  v'.  Point  S  is  the  zero 
intercept  of  Line  v.  The  iteration  is  stopped  here  because  the  difference  between  the 
Points  R  and  S  is  sufficiently  small. 

is  determined  by  identifying  the  value  of  which  gives  Al7p=0.  The  initial 
value  used  to  determine  is  found  in  a  manner  similar  to  that  for  given  by  Eq. 
(21).  Because  the  sign  of  q  is  changed,  there  are  sign  changes  in  Eq.  (21),  giving 
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Figure  9  Slope  of  Curve  B  in  Figure  8.  The  bold  line  is  the  actual  slope,  and  the  thin  solid  lines 
are  approximations  using  Eq.  (IS).  The  dashed  lines  exaggerate  the  error  between  the 
approximations  and  actual  slope  by  a  factor  of  500. 


Algorithm  1 1  shows  the  implementation  of  the  iterative  search  used  to  find  for  a 
position  given  as  y-y^,  L~Ly,  and  sin(<}>)—St,.  If  Lu<Lo  then  the  location,  L^,  s„,  is 
outside  of  Xq  if  U  for  ions  is  greater  than  U  of  Xq.  This  can  be  seen  in  Figure 
1  where  the  contour  lines  are  equal  to  U  since  q  is  positive.  If  the  measured  position  is 
outside  of  Xq,  then  is  given  the  value  of  -1  to  indicate  that  is  not  defined.  The 
initial  upper  value  of  for  the  search  is  given  by  Eq.  (24).  L<L^  because  is  only 
defined  for  positions  inside  of  Xq.  Substitution  of  Eq.  (7)  into  Eq.  (24)  shows  that  Eq. 
(24)  gives  a  positive  result  when  LKLq.  The  initial  lower  value  is  set  to  zero.  A  line 
is  used  to  approximate  the  relationship  between  AUp  and  using  two  points  and  is  used 
to  solve  for  a  new  value  of  The  point  farthest  from  AUp=0  is  replaced  by  the  new 
solution  for  for  the  next  iteration.  The  iteration  is  continued  until  the  value  of 
W^=  Wjy  is  known  to  sufficient  accuracy,  usually  in  about  three  iterations. 

Figure  10  illustrates  the  solution  for  for  the  conditions  given  for  Figures  7 
and  8,  and  Points  A  through  D  have  the  same  energy  as  do  those  in  the  previous  figures. 
The  X  axis  shows  the  y  axis  shows  AUp,  and  the  dark  curve  shows  the  relationship 
between  the  two.  The  curve  is  shown  as  a  dotted  line  on  the  right  side  of  the  figure 
where  a  peak  does  not  exist  at  AUp  such  as  for  that  shown  by  Curve  D  in  Figure  8.  The 
initial  upper  boundary  at  Point  C  is  identified  using  Eq.  (24),  and  the  initial  lower 
boundary  at  Point  A  is  at  W„=0.  The  thin  solid  line  connecting  Points  A  and  C  shows 
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Algorithm  11  Calculation  of  ion  lower  cutoff  energy  for  a  position  givoi  by  L,  y,  and 
sin(^). 

FUNCTION  DEFINITION  OF  W^ly,  L,  s) 
local  constants  required  for  W|.(y,  L,  s): 
k  HE  Bo'wa^ 

local  variables  required  for  W^ly,  L,  s): 

Uo,  J„  W.  Wo.  W„  AU.  AUo.  AU, 
implementation  of  W^ly.  L,  s): 

Uy  =  Ar(y.  1 .  U  Value  of  li  at  y,L,s  if  W=  1  eV;  p  =  PrW  Egn  3 

J,  =  J(y,  1,  L)  Value  of  J  at  y,L,s  if  W-  1  eV;  J  =  Eqn  4 

Uo  =  Ulfalse,  0,  L,  s)  Value  of  U  at  y,L,s  if  W=  0  eV;  U  =  Uo+W  Egn  5 

Wo  =  0 

AUo  =  AUlO.  0,  Uo,  false,  U,  -1) 

if  (L  >  Lo)  or  [AUo  >  01  then  satellite  is  outside  of  3£o 

return  value  of  -1  and  exit  function  -1  indicates  that  does  not  exist 
W,  =  W,.^-(y  a-A-L^'  -  knJ)  /  aW/aL(/y,-W„„.  L)  Equation  24 

AU,  =  AUOi/i-W,,  Ji-W,*,  Uo  +  W,,  false.  PeakWil/;,  W„  Ji-W,*),  -1)  * 
while  abs|W,-Wo)  >  +  repeat  following  loop: 

begin  while  hop 

VI  =  Wo-  {W,-Wo)/|AU,-AUo)-AU 

AU  =  AUU/,-W,  J,-W»‘,  Uo  +  W,  false,  Peak.W^_{p,^W,  J,-W*),  -11  * 
if  abs(AU,)  <  abslAUol  then 
begin  conditional  block 
Wo  =  W 
AUo  =  AU 
end  conditional  block 

else 

begin  conditional  block 
W,  =  W 
AU,  =  AU 
end  conditional  block 
end  while  hop 
return  value  of  W 

*  These  lines  can  be  made  more  efficient  by  calculating  p,-W„  and  only  once 

for  each  line  instead  of  twice. 

the  linear  approximation  to  the  relationship.  The  zero  intercept  of  this  approximation  is 
very  close  the  actual  solution,  shown  by  Point  B.  The  process  converges  to  Point  B  in 
one  to  two  more  iterations. 


6.  CALCULATION  OF  ION  UPPER  CUTOFF  EP^RGY, 

Figure  1 1  illustrates  the  ion  upper  cutoff  energy,  W„.  Six  sets  of  ion  drift  paths 
in  the  equatorial  plane  are  shown  in  this  figure,  and  these  are  calculated  using  the  method 
presented  in  Section  3  and  Algorithm  6.  The  constants  of  motion  are  determined  for  ions 
with  one  of  six  different  values  of  measured  at  the  white  circle  with  iv=3.5  Re, 
y.„=0.90,  and  sin(<^.„)=0.95.  In  each  case,  W^.  The  six  values  of  are  given 
on  the  right  side  of  the  figure  and  are  labeled  A  to  G.  The  letter  C  is  omitted  so  that  the 
energies  shown  in  this  figure  will  match  the  energies  used  in  following  figures.  In  this 
example  Io=6  Rg,  and  the  shaded  region  is  bounded  by  Line  A  shows  the  path  for 
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Figure  10  Relationship  between  AU,  and  for  the  conditions  given  for  Figure  7.  This 
illustrates  the  solution  for 

ions  with  Wj^=5.4  keV.  In  this  case  there  is  a  closed  westward  path  around  the  Earth 
inside  of  Xq,  and  there  is  a  separate  sunward  path  open  to  the  tail  on  the  dawn  side  of 
the  Earth.  This  path  is  similar  to  Lines  C  and  D  of  Figure  7.  The  point  of  measurement 
is  on  the  path  open  to  the  tail.  Line  B  shows  the  path  for  ions  with  W^=9.4  keV.  This 
energy  is  labeled  Wg  to  indicate  the  stagnation  energy  which  is  discussed  in  Section  7. 
Here  the  iimer  closed  westward  path  joins  with  the  outer  open  path  at  a  single  point  on 
the  dawn  meridian  near  3.2  Re,  and  the  point  of  measurement  is  located  on  the  outer 
open  path.  Line  D  shows  the  path  for  ions  with  W^=10.8  keV.  Here  the  inner 
westward  path  is  joined  with  the  outer  open  path  to  form  a  single  open  path  that  crosses 
the  dawn-dusk  meridian  on  the  dusk  side  of  the  Earth.  Line  E  shows  the  path  for  ions 
with  W„.  This  path  is  similar  to  Line  B  except  that  the  point  of  measurement  is 
now  located  on  the  iimer  closed  path.  This  is  the  transition  trajectory  between  open  and 
closed  drift  paths,  because  all  ions  with  are  on  open  paths  while  higher 

energies,  as  shown  by  Lines  F  and  G,  place  the  point  of  measurement  on  closed 
westward  paths.  In  the  last  two  cases,  the  open  paths  from  the  tail  are  separate  from  the 
open  sunward  paths. 

Wg  is  calculated  in  a  manner  similar  to  except  that  AU  is  evaluated  to  find 
a  single  solution  at  dawn  instead  of  dusk.  Figure  12  shows  AU  as  a  function  of  L-shell 
at  dawn  for  the  conditions  given  for  Figure  1 1 .  Algorithm  6  is  used  to  calculate  AU,  and 
each  curve  shows  AU  versus  L-shell  at  dawn  for  the  same  six  values  of  W*,  and  for  the 
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Figure  11  Convection  paths  of  singly  diarged  ions  measured  at  the  white  circle  for  six  values  of 
W^f.  The  shaded  region  is  bounded  by  and  arrows  show  the  direction  of  drift  along  the  paths. 


same  measured  location  given  for  Figure  11.  A  seventh  energy  is  included  for  Curve  C. 
The  peaks  in  this  figure  look  very  similar  to  those  seen  in  Figure  4  for  electrons  at  dusk. 
Curve  E  shows  the  solution  for  W„  because  it  gives  a  single  solution  of  AU=0  at  dawn 
while  placing  the  point  of  measurement  on  the  inner  closed  path.  The  peak  values  of  AV 
at  dawn  are  referred  to  as  AUp. 

In  order  to  determine  the  value  of  AUp,  the  L-shell  Lp  at  dawn  at  which  the  peak 
occurs  is  identified,  so  that  AUp—AU(Lp).  Lp  occurs  at  the  L-shell  where  dAUIdL—O  at 
dawn.  AU=Um-U,  and  U  comes  from  Eq.  (5)  with  q=  1  and  1,  so  that 

^(L=LJ  =  -.^(L=L.)  -  yuALr'  =  0  (25) 

aZ,  ^  dL  ^  ® 


This  solution  is  identical  to  Eq.  (22)  exc^t  for  the  sign  change  in  the  second  term.  It 
is  solved  iteratively  as  is  Eq.  (11)  by  using  Eq.  (15),  exc^t  that  the  solution  of  interest 


Figure  12  MJ  as  a  function  of  L-shell  at  dawn  (sin(<^)=l)  for  the  same  conditions,  measured 
location,  and  energies  as  those  given  for  Figure  11. 


is  that  given  by  Eq.  (17)  where  Q^=^aA,  and  Qc--C2-  The  values  of 

c,  and  C2  of  Eq.  (IS)  are  first  determined  from  initial  points  at  and  L=Io-l.  An 
estimate  of  Lp  is  determined  using  Eq.  (17),  and  this  new  value  is  used  to  replace  the 
initial  estimated  value  of  L-shell  farthest  from  the  new  estimate.  New  values  for  Cj  and 
C2  are  determined,  which  then  give  a  better  estimate  of  Lp.  This  process  is  repeated  until 
it  converges  to  a  solution,  usually  in  four  to  five  iterations.  Algorithm  12  shows  this 
method  of  determining  Lp. 

Figure  13  illustrates  Algorithm  12  for  the  solution  of  Lp  for  the  case 
Curve  E  of  Figure  12.  The  bold  line  shows  dAUIdL  as  a  function  of  L-shell  at  dawn. 
Lp  is  located  where  this  line  crosses  zero  (Point  U).  Points  P  and  Q  at  L=Lo  and  L—L^-X 
are  used  to  determine  initial  values  for  c,  and  Cj  [Eq.  (IS)].  Line  u,  which  is  a  thin  solid 
line,  shows  the  approximation  to  dAUIdL  by  using  Eq.  (IS).  Because  the  approximation 
is  close  to  the  actual  solution.  Line  u',  which  is  a  thin  dashed  line,  is  shown  with  the 
error  between  the  approximation  and  the  actual  solution  exaggerated  by  a  factor  of  SO. 
Equation  (17)  is  used  to  identify  the  zero  intersect  of  Line  u  at  Point  R.  The  second 
iteration  uses  Points  Q  and  R  to  give  Line  v  which  yields  Point  S.  This  method  is 
iterated  two  more  times  to  yield  Points  T  and  U  before  it  stops. 

Wg  is  determined  by  identifying  the  value  of  which  gives  AUp—Q.  The  initial 
value  used  to  determine  W„  is  found  in  a  manner  similar  to  that  for  W^.  given  by  Eq. 
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Algorithm  12  Calculation  of  Lf  at  dawn  for  given  values  of  fi  and  J  for  ions. 


FUNCTION  DEFINITION  OF  PeakW„(;/,  J) 
local  constants  required  for  PeakW^Ou,  J); 


k  »  -Bo-a^-o; 


local  variables  required  for  PeakW^Ov,  J): 

L,  L,,  Lj  These  hold  the  current  solutions  of  L-shetl 

D,  D,,  Dj  These  hold  dW/dL  for  the  current  solutions  of  L-shell 

Cf,  Cj,  Qg,  Qq 
implementation  of  PeakWH(/y,  J): 

Li  =  Lo 

Di  =  aW/3L|//,  J,  L,) 

L2  =  L,  -  1.0 

Dj  =  aw/aiu/,  J,  Lj) 

Q*  =  -K-a-A 

while  abslLj-L,)  >  6^  repeat  the  following  loop: 
begin  while  loop 
Cj  =  (L,^-D,  - 

Cl  =  -  Cj'L, 

Qg  =  k  -  Cl 


Find  two  starting  points  for  the  iteration. 

L0  is  usually  within  a  few  of  the  final  solution. 


First  coefficient  for  Equation  1 7. 


■  llj-^’) 


Qq  —  Qft  +  4-‘Q^'C2 


L  =  [%(-Qb  - 
D  =  3W/aL|//,  J,  L) 
if  abs(L  -  L,)  <  abs(L  -  Lj)  then 
begin  conditional  block 
L,  =  L 
Dj  =  D 

end  conditional  block 

else 

begin  conditional  block 
L,  =  L 
D,  =  D 

end  conditional  block 
end  while  loop 
return  value  of  L 


Solution  to  Equation  15 
coefficients. 

Second  coefficient  for  Equation  1 7. 
Value  inside  of  root  in  Equation  1 7. 
Equation  1 7  yields  better  value  for  L. 
New  value  for  dW/dL  is  computed. 

Replace  whichever  solution  is 
farthest  from  new  solution. 


(21).  Because  the  sign  of  q  is  changed  and  the  sign  of  sin(<^)  is  changed,  there  are  some 
sign  changes  in  Eq.  (21),  giving 


Algorithm  13  shows  the  implementation  of  the  iterative  search  used  to  find  for  a 
position  given  as  y=yu,  and  sin(«^)=Sj,^.  The  initial  lower  value  of  for  the 

search  is  given  by  Eq.  (26).  TTiis  equation  always  gives  a  positive  solution.  TTie  initial 
upper  value  is  chosen  arbitrarily  to  be  10  keV  greater  than  the  lower  value.  A  third  point 
is  picked  by  assuming  that  the  first  two  points  lie  on  parabola,  and  that  the  slope  of  the 
parabola  is  zero  at  the  lower  point.  A  parabola  is  used  to  approximate  the  relationship 
between  bJJp  and  W^,  using  three  points,  and  it  is  used  to  solve  for  a  new  value  of 
The  point  that  is  farthest  from  AUp=0  is  replaced  by  the  new  solution  for  for  the 
next  iteration.  The  iteration  is  continued  until  the  value  of  is  known  to 
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Figure  13  Slope  of  Curve  E  in  Figure  12.  The  bold  line  is  the  actual  slope,  and  the  thin  solid 
lines  are  approximations  using  Equation  IS.  The  dashed  lines  exaggerate  Ae  error  between  the 
approximations  and  actual  slope  by  a  factor  of  50. 


sufficient  accuracy,  usually  in  about  three  iterations. 

Figure  14  illustrates  the  solution  for  for  the  conditions  given  for  Figure  13. 
Points  A  through  G  have  the  same  energy  as  those  in  the  previous  figure.  The  x  axis 
shows  W^,,  the  y  axis  shows  AUp,  and  the  dark  curve  shows  the  relationship  between  the 
two.  The  initid  lower  boundary  at  Point  D  is  identified  using  Eq.  (26),  and  the  initial 
upper  boundary  at  Point  G  is  found  by  adding  10  keV  to  the  energy  of  Point  D.  Point 
F  is  found  using  the  zero  intercept  of  the  parabola  formed  from  Points  D  and  G  with  a 
slope  of  zero  at  Point  D  (shown  by  thin  dotted  line).  The  zero  intercept  of  the  parabola 
formed  from  these  three  points  (thin  solid  line)  is  very  close  to  the  true  solution  shown 
by  Point  E.  The  process  converges  to  Point  E  on  the  next  iteration. 


7.  CALCULATION  OF  ION  STAGNATION  ENERGY,  Wg 

The  ion  stagnation  energy,  Wj,  is  the  energy  for  which  ions  on  open  drift  paths 
switch  from  crossing  the  dawn-dusk  meridian  on  the  dawn  side  to  crossing  on  the  dusk 
side.  This  energy  is  called  a  stagnation  energy  because  ions  with  this  energy  move 
through  a  stagnation  point  on  the  dawn  side  of  the  Earth,  even  though  they  may  still  be 
on  open  drift  paths.  The  significance  of  Wg  is  discussed  in  Kerns  et  al.  [1993].  Figure 
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Algorithm  13  Calculation  of  ion  upper  cutoff  energy,  Wh,  for  a  position  given  by  L,  y,  and 
sin(^). 


FUNCTION  DEFINITION  OF  W„|y,  L,  s) 
local  constants  required  for  WmIv,  L,  s): 
k  *  Bo-wa^ 

local  variables  required  for  WhIy,  L.  s): 

Uo,  Hu  J„  Wo,  W„  Wj,  AUo,  AU„  AUj,  Q^,  Ob,  Qc 
implementation  of  Wh(v,  L,  s): 

H\  -  HiV>  1 .  U  Value  of  fj  at  yA,s  if  W=  1  eV;  ij  =  ;/,-W  Eqn  3 

J,  =  Jly,  1 ,  U  Value  of  J  at  y,L.s  if  W=  1  eV;  J  =  J,  W^  Eqn  4 

Uo  =  U(false,  0,  L,  s)  Value  of  U  at  y.L.s  if  W-0  eV;  U  =  Uo+  W  Eqn  5 
Wo  =  +  K-a-A-L*^’)  /aW/aL0i,-W„^,  J,-W^’‘,  L)  Equation  26 

AUo  =  AUUy,-Wo,  J^Wo^  Uo  +  Wo,  false.  PeakWnU^rWo,  JrWo’‘l,  1}  * 

The  second  starting  point  is  offset  from  the  first  by  W„,^. 

W,  =  Wo  +  W„,p 

AU,  =  AU(iU,-W,.  J,  W/,  Uo  +  W,.  false.  PeakW„0L/,  W„  J,  W,’‘),  1)  * 

A  parabolic  approximation  is  used  to  find  the  third  point.  It  uses  the  points 
and  W,j  and  assumes  the  slope  at  Wg  is  zero. 

Wj  =  Wo  +  W„,p-[AUo/(AUo-AU,)]’‘  This  is  zero  intercept  of  parabola. 

AUj  =  AUlPi-Wj.  Ji-Wj'^,  Uo  +  Wj.  false.  PeakWH(;/,-W2,  JrWj**),  1)  * 
Reorder  Wg,  W,.  and  W,  so  that  \b-Ug\  <  |  At/,|  <  lAt/^l 
if  abs(AUi)  <  abs(AUo)  then 

swap  values  of  AU,  with  AUo  and  W,  with  Wo 
if  abslAUj)  <  abslAU,)  then 

swap  values  of  AUj  with  AU,  and  Wj  with  W, 
if  abs(AU,)  <  abs(AUo)  then 

swap  values  of  AU,  with  AUo  and  W,  with  Wo 
while  abs<W,-Wo)  >  tfw'IWo  +  W,)  repeat  following  loop: 
begin  while  loop 

0^  =  (AUo-|Wj-W,)  +  AU,-|Wo*W2)  +  AU2-{W,-Wo)]  / 

(Wo*-(Wj-W,)  +  W,2-|Wo-W2)  +  W22-(W,.Wo)] 

Qa  =  [AUo-AU,  -  Q*-(Wo*-W,")]  /  (Wo-W,) 

Oc  =  AUo  -  Wo-IQb  +  Q^  Wo) 

AUz  =  AU,  Shift  point  1  to  point  2 

Wj  =  W,  and  point  0  to  point  1, 

AU,  =  AUo  and  find  new  point  0. 

W,  =  Wo 

Wo  =  ’/aT-Qs  +  (Qb^  -  4-Q^-Qc)*]  /  Qa  Quadratic  Formula 

AUo  =  -ttU{;/,  Wo.  J,-Wo’‘,  Uo  +  Wo,  false,  PeakW„(^,-Wo,  J,-Wo’‘>.1)  * 
end  while  loop 
return  value  of  Wo 

*  These  lines  can  be  made  more  efficient  by  calculating  p,-W„  and  J,'W„^  only  once 
for  each  line  instead  of  twice. 


1 1  illustrates  the  ion  stagnation  energy,  In  this  figure  Lines  A,  B,  and  D  place  the 
point  of  measurement  on  an  open  drift  path.  Line  A,  which  has  W,|^<  Wj  crosses  the 
dawn-dusk  meridian  on  the  dawn  side  of  the  Earth.  Line  D,  which  has  Wj  crosses 
the  dawn-dusk  meridian  on  the  dusk  side  of  the  Earth.  Line  B  shows  Vi^=  Wj,  and  here 
the  closed  inner  path  meets  the  open  outer  path  at  a  single  point  at  dawn.  Although  the 
point  of  measurement  is  on  an  open  field  line,  ions  convecting  from  the  tail  move 
through  the  dawn  stagnation  point  to  reach  the  point  of  measurement,  and  cannot  convect 
to  the  point  of  measurement. 
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W^.  keV 

Figure  14  Relationship  between  AUp  and  for  the  conditions  given  for  Figure  13.  This 
illustrates  the  solutions  for  and 

Wg  is  calculated  in  a  manner  identical  to  that  of  Wg  except  that  Eq.  (26)  gives  the 
initial  upper  value  instead  of  the  initial  lower  value  of  V/y  in  the  search.  Also,  the 
quadratic  formula  uses  the  negative  solution  to  find  &e  zero  intercept  of  the 
approximating  parabolas.  The  upper  solution  to  AXJp=0  is  shown  by  Curve  E  in  Figure 
12  where  Wy=  Wg.  The  lower  solution  to  AUp=0  is  shown  by  Curve  B,  and  here 
Wy^  Wg.  Curve  D  uses  the  value  given  by  Eq.  (26)  for  Wy,  and  as  can  be  seen  in  the 
figure,  Wj  is  less  than  this  value.  The  location  of  the  peaks  are  determined  using 
Algorithm  12. 

Algorithm  14  shows  the  implementation  of  the  iterative  search  used  to  find  Wg 
for  a  position  given  as  y=yy,  L=Ly,  and  sin(<^)=Sjj,.  The  initial  upper  value  of  Wy  for 
the  search  is  given  by  Eq.  (26).  This  equation  always  gives  a  positive  solution.  The 
initial  lower  value  is  set  to  one-half  the  initial  upper  value.  A  third  point  is  picked  by 
assuming  the  first  two  points  lie  on  parabola,  and  the  slope  of  the  parabola  is  zero  at  the 
upper  point.  If  Wy  for  the  third  point  is  less  than  or  equal  to  zero,  then  Wy  for  this  point 
is  set  to  one-half  the  lower  value.  A  parabola  is  us^  to  approximate  the  relationship 
between  AUp  and  Wy  using  three  points,  and  it  is  used  to  solve  for  a  new  value  of  Wy. 
The  point  which  is  farthest  from  AUp=Q  is  r^laced  by  the  new  solution  for  Wy  for  the 
next  iteration.  The  iteration  is  continued  until  the  value  of  Wg=  Wy  is  known  to  sufficient 
accuracy,  usually  in  about  three  iterations. 
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Algorithm  14  Calculation  of  ion  stagnation  energy,  Wj,  for  a  position  given  by  L,  y,  and  8in(^). 


FUNCTION  DEFINITION  OF  W,(y,  L,  s) 
local  constants  required  for  W^ly,  L,  s): 
k  s  Bo'wa* 

local  variables  required  for  W,(y.  L,  s); 

Uo.  a,.  Jv  Wo.  W„  Wj,  AUo,  AU„  AUj,  Q^,  Qe,  Qc 
implementation  of  Wg(y,  L,  s): 

//,  =  //(y,  1 ,  L)  Value  of  fj  at  y,L,s  if  W=  1  eV;  =  p,-W  Eqn  3 

J,  =  J(y,  1 .  U  Value  of  J  at  y,Ls  if  1  eV;  J  --  J,-W^  Eqn  4 

Uo  =  U(false.  0,  L,  s)  Value  of  U at  y,Ls  if  W=0  eV;  U  =  Uo+W  Eqn  5 
Wo  =  -W„^  (k/L^  +  K-a  A-L'^n  /aW/aL(;/,  W„^,  J,•W«^^  L)  Equation  26 
AUo  =  AU0;,-Wo,  Ji  Wo*^,  Uo  +  W'o,  false,  PeakW„|/;,  Wo.  J^Wo*),  1)  ★ 

The  second  starting  point  is  offset  from  the  first  by  W^.  - 
W,  =  ’AWo 

AU,  =  AU(/;i  W„  J,•W,^  Uo  +  W,,  false,  PeakW„(iU,  W„  J,  W,’‘),  1)  * 

A  parabolic  approximation  is  used  to  find  the  third  point.  It  uses  the  points  Wg 
and  W,j  and  assumes  the  slope  at  Wg  is  zero. 

Wj  =  Wo  -  (Wo-W,)-(AUo/(AUo-AU,)]’‘  This  is  zero  intercept  of  parabola. 

if  W2  s  0  then 

Wj  =  ’AW, 

AUj  =  AU0/,  W2.  Ji-Wj'*,  Uo  +  Wj.  false,  PeakW„0;,-W2,  Ji  Wj’^),  1)  ★ 
Reorder  VJg,  W,,  and  VJ^so  that  |Ai/o|  <  lAty,!  <  |At41 
if  abs(AU,)  <  abslAUo)  then 

swap  values  of  AU,  with  AUo  s^d  W,  with  Wq 
if  abslAUj)  <  absiAU,)  then 

swap  values  of  AUj  with  AU,  and  W^  with  W, 
if  abs{AU,)  <  abslAUo)  then 

swap  values  of  AU,  with  AUo  and  W,  with  Wo 
while  abs|W,-Wo)  >  dw'iW^o  +  W,)  repeat  following  loop: 
begin  while  loop 

Qa  =  [AUolWj-W,)  +  AU,-{Wo-W2)  +  AU2(W,-Wo)1  / 

[Wo'-lWj-W,)  +  W,2-(Wo-W2)  +  W22.(W,-Wo)] 

Qb  =  [AUo-AU,  -  Qa-IWo^-W,^)]  /  IWo-W,) 

Qc  =  AUo-Wo(Qb  +  Qa-Wo) 

AUj  =  AU,  Shift  point  1  to  point  2 

Wj  =  W,  and  point  0  to  point  1, 

AU,  =  AUo  and  find  new  point  0. 

W,  =  Wo 

Wo  =  ’A  l-Qe  *  ' 4  Qa  Qc)’‘1  /  Qa  Quadratic  Formula 

if  Wo  s  0  then 

Wo  -  %W, 

AUo  =  AUU/,’Wo.  J,  Wo*.  Uo  +  Wo,  false,  PeakW„{A/,  Wo.  J,  Wo’‘),1)  * 
end  while  loop 
return  value  of  Wo 

*  These  lines  can  be  made  more  efficient  by  calculating  p,-W„  and  only  once 

for  each  line  instead  of  twice. _ 

Figure  14  illustrates  the  solution  for  Wg  for  the  conditions  given  for  Figure  13. 
The  initial  upper  boundary  at  Point  D  is  identified  using  Eq.  (26),  and  the  initial  lower 
boundary  at  Point  A  is  one-half  the  energy  of  Point  D.  Point  C  is  found  using  the  zero 
intercept  of  the  parabola  formed  from  Points  D  and  A  with  a  slope  of  zero  at  Point  D 
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(shown  by  thin  dotted  line).  The  zero  intercept  of  the  parabola  formed  from  these  three 
points  (thin  solid  line)  is  very  close  to  the  U\ie  solution  shown  by  Point  B.  The  process 
converges  to  Point  B  on  the  next  iteration. 
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